#ifndef TREEFORCEFINDER_H
#define TREEFORCEFINDER_H

#include <iostream>
#include <stdio.h>
#include <deque>
#include <list>

#include "MACRO.H"
#include "Vector.h"
#include "Particle.h"
#include "Units.h"

	
class TreeForceFinder {	
	public:
	
		double theta;
		double eps;
		
		Units u;
		
		TreeForceFinder(): 
		theta(0.6), 
		eps(0.02) 
		{ 
			cell0.S  = Vector<double>(100,100,100);
			cell0.R0 = Vector<double>(-50,-50,-50);
		}
		
		class Cell {
			public:
	
			static int n;	
				
			double m; // mass
			
			Vector<double> R0; // left bottom corner
			Vector<double> S;  // size
			Vector<double> CM; // center of mass
	
			vector<Cell> cell;
	
			vector< Particle* > p;
			
			Cell() : m(0) {}
		
			void open(){ 

				Vector<double> Ss = S*.5;
				Vector<double> C = R0 + Ss;
					
				{
					
				vector<Cell> cell(8);
				
				for (int i=0; i<p.size(); i++)
					cell[	(p[i]->R.x < C.x ? 0 : 1) + 
							(p[i]->R.y < C.y ? 0 : 2) + 
							(p[i]->R.z < C.z ? 0 : 4)].p.push_back(p[i]);

				p.clear(); vector<Particle*>().swap(p);		

				
				for (int i=0; i<2; i++)
				for (int j=0; j<2; j++)
				for (int k=0; k<2; k++)
					cell[i+j*2+k*4].R0 = Vector<double>(	i==0 ? R0.x : C.x, 
																		j==0 ? R0.y : C.y, 
																		k==0 ? R0.z : C.z);
				
				int size(0);
				for (int i=0; i<cell.size(); i++)
				if (!cell[i].p.empty()) size++;
				
				this->cell.resize(size); size=0;
				for (int i=0; i<cell.size(); i++)
				if (!cell[i].p.empty())
					this->cell[size++] = cell[i];
				
				}

				
				for (int i=0; i<cell.size(); i++) {
					
					cell[i].S = Ss;
				
					Vector<double> CMs;
					double ms(0);
					for (int j=0; j<cell[i].p.size(); j++) {
						CMs+= cell[i].p[j]->m*cell[i].p[j]->R;
						ms += cell[i].p[j]->m;
					}
					cell[i].CM = CMs/ms;
					cell[i].m = ms;
					
					if (1<cell[i].p.size()) cell[i].open();
						else n++;
							
					vector<Particle*>().swap(cell[i].p);
				}
				
			}
			
		
			void calc_acc(Vector<double> &R, Vector<double> &A, double &theta, double &eps) {
								
				Vector<double> D = R-CM;
				double d2 = D.sqr() + SQR(eps);
		
				if (d2==0.0) return;
			
				double d = sqrt(d2);
				double angle = S.x/d;				
								
				if ( theta<angle && cell.size() )
				for (int i=0; i<cell.size(); i++) 
					cell[i].calc_acc(R, A, theta, eps);
				else {
					double mZd3(m/(d*d2));
					A-=D*mZd3;
				}
				
			}
		
			void print() {
				
				printf("\n");
				printf(_B"TreeForceFinder::cell[%u]\n"B_, (long unsigned int*)this);
				printf(" ├ m = %g\n", m);
				printf(" ├ R0 = {% 12.8e, % 12.8e, % 12.8e}\n", R0.x, R0.y, R0.z );
				printf(" ├ S  = {% 12.8e, % 12.8e, % 12.8e}\n", S.x, S.y, S.z );
				printf(" └ cell.size() = %i\n", cell.size());
				/*
				printf(" cell = {");
				for (int i=0; i<cell.size(); i++)
				printf("cell[%u]%s", (long unsigned int*)(&cell[i]), i<cell.size()-1 ? ", " : "");
				printf("}\n");
				*/
				printf("\n");
		
			}
			
			
		
	} cell0;
	
	
	void calc_accs(vector<Particle> &p) {
		
		for (int i=0; i<p.size(); i++) {
			p[i].A*=0;
			calc_acc(p[i].R, p[i].A);
		}
		
	}
	
	void calc_acc(Vector<double> &R, Vector<double> &A) {
		
		cell0.calc_acc(R, A, theta, eps);
		
		A*=u.G();
	}
	
	void build_tree(vector<Particle> &p){
		
		cell0.m = 0;
		vector<Cell>().swap(cell0.cell);
		double rmax(0); 
		
		for (int i=0; i<p.size(); i++) {
					
			if ( SQR(rmax) < p[i].R.sqr() ) 
				rmax = p[i].R.norm();
			
			cell0.m += p[i].m;
			cell0.p.push_back(&p[i]);
			
		}
		
		cell0.open();
		
	}

	void test_tree(vector<Particle> &p) {
		
		Vector<double> R(20,0,0);
		Vector<double> A; 

		calc_acc(R, A);
		
		A.print();
		
		Vector<double> A2; 
		for (int i=0; i<p.size(); i++) {
			
				Vector<double> DR = R - p[i].R;
				double dr2 = DR.sqr();
				double dr = sqrt(dr2);
				
				A2-= (p[i].m/(dr2*dr))*DR;
				
		}
		A2*=u.G();
		
		A2.print();
		
	}
	
	void free(){
		printf("TreeForceFinder::free()...\n");
		vector<Cell>().swap(cell0.cell);
	}
	
	
	void print(){
		printf("TreeForceFinder::print()\n");
		printf("  Cell::n = %i\n", cell0.n);
	}
	
	
};
	
int TreeForceFinder::Cell::n = 0;
	
	
#endif















